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Abstract 


In  this  report  we  present  test  results  on  various  non-conventional  signal  processing  approaches 
as  applied  to  data  from  regional  seismic  arrays.  In  particular,  an  adaptive  approach  for 
automatic  identification  of  regional  arrivals  was  tested  to  ( Der  et  al  1993).  This  method  was 
designed  to  circumvent  some  difficulties  that  plagued  previous  attempts  of  3-component 
polarization  analyses  for  regional  arrival  identification.  In  particular,  the  methodology  is 
designed  to  avoid  complications  due  to  signal  distortion  by  near  sensor  geology  and  complex 
arrivals  made  up  of  superpositions  of  multiple  independent  signal  components.  The  idea  is  to 
replace  time  delay  operators  and  assumed  particle  motion  patterns  used  in  deterministic  signal 
processing  with  empirically  computed  sets  of  transfer  functions  using  a  learning  set  of  relatively 
high  S/N  events.  The  application  of  this  approach  requires  more  than  three  sensor  outputs  (such 
as  a  3-component  combination  plus  some  extra  array  sensors). 

We  have  tested  both  a  frequency  domain  and  a  time  domain  approach  to  this  algorithm.  Thus  far 
the  frequency-domain  version  has  reached  a  further  stage  of  development.  The  methodology 
appears  to  be  an  attractive  alternative  to  conventional  array  processing  methods  in  quickly 
identifying  and  discriminating  various  arrivals  from  the  same  source  region,  or  similar  arrivals 
from  different  source  regions  at  a  small  array. 

We  are  also  addressing  the  problem  of  onset  time  estimation  for  regional  arrivals.  Most  methods 
proposed  thus  far  attempt  to  find  the  first  point  where  the  statistical  properties  of  the  seismic 
trace  change.  This  is  difficult  or  next  to  impossible  with  the  often  emergent  regional  arrivals. 
Instead,  we  are  advocating  a  methodology  based  on  the  cumulative  sum  (CUSUM)  of  various 
statistics.  Superposing  a  linear  trend  on  the  cumulative  sum  can  make  the  changes  in  the  trend 
of  CUSUM  to  become  minima  of  the  resulting  function.  This  function  lends  itself  to  the 
application  of  numerous  algorithms  for  finding  the  minima  of  a  function.  Besides  the  viable 
method  of  onset  time  estimation  by  finding  of  the  absolute  minima,  repeated  applications  of 
simulated  annealing  (SA)  results  in  populations  of  onset  time  estimates  that  can  be  treated  with 
ordinary  statistical  methods.  Tests  of  this  methodology  are  very  promising.  CUSUM-based 
onset  time  estimation  for  Pn  appears  to  be  competitive  to  estimation  of  onset  times  by  an 
analyst. 

Key  Words:  seismology,  arrays,  signal  processing,  statistics,  detection,  identification,  travel 
time. 


IV 


Objective 


The  objective  of  the  work  thus  far  is  the  facilitation  of  the  automatic  recognition  of  various 
regional  arrivals.  This  is  accomplished  by  a  computational  scheme  for  the  adaptive  learning  of 
the  generalized  particle  motion  patterns  as  seen  by  small  arrays  that  may  also  contain  three- 
component  sensors. 

The  purpose  of  the  work  is  to  relieve  the  human  operator  in  seismic  monitoring  from  the  effort 
to  identifying  and  locating  routine  blasting  at  various  quarries  and  to  enable  him  to  possibly 
identify  the  sources  of  seismic  waves  from  the  Lg  phase  alone  and  to  point  out  unusual  events 
near  some  quarries  that  may  merit  further  investigation.  The  techniques  proposed  do  not  have  to 
be  applied  alone,  but  can  be  used  jointly  with  other  existing  methods. 

In  the  later  part  of  this  report  we  present  results  obtained  by  applying  a  new  method  that 
combines  cumulative  sums  of  a  test  statistic  (CUSUM)  and  simulated  annealing  (SA),  to  solve 
the  problem  of  automatic  onset  time  estimation  for  regional  seismic  arrivals.  Onset  time 
estimation  and  segmentation  of  a  regional  seismogram  into  various  kinds  of  seismic  phase  wave 
groups  is  part  of  an  analyst’s  routine  which  appears  to  be  amenable  to  automation. 

Organization  of  this  Report 

This  report  consists  of  three  parts.  In  the  first  part,  we  describe  the  testing  of  generalized 
polarization  processor  (Der  et  al  1993).  In  the  second,  we  discuss  the  development  of  an  onset 
time  estimator  based  on  CUSUM  and  SA.  Finally,  in  the  third  part  we  describe  work  on  other 
topics,  in  particular,  modification  made  to  the  deconvolution  program. 
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ANALYSES  OF  THE  GENERALIZED  POLARIZATION  PROCESSOR 


Introduction 

Identification  of  regional  “phases”  at  small  seismic  arrays,  as  currently  practiced,  is  largely 
based  on  their  slownesses  and  rectilinearity,  frequency  content  and  the  time  context  they  appear 
in  on  the  seismograms.  The  concepts  of  slowness  and  rectilinearity  implicitly  imply  the 
assumptions  of  plane  waves  and  linear  particle  motion  as  models  of  seismic  signals. 
Rectilinearity  is  based  on  simple  halfspace  models  of  P  and  S  motion  that  have  little  relevance 
to  the  real  world.  Although  these  models  fit  many  situations  reasonably  well,  there  are  too  many 
cases  where  such  simple  models  fit  the  data  only  poorly.  The  issue  is  the  quality  of  these  fits. 
The  question  arises  whether  one  could  construct  some  more  complex  signal  models  that 
describe  the  spatial  variations  of  waveforms  better.  Moreover,  if  such  models  can  be  found  they 
may  have  onsiderable  discriminative  value  in  distinguishing  the  various  arrival  types  at  regional 
arrays.  Assuming  more  complex  multicomponent  structure  for  regional  waveforms  has  been 
shown  to  fit  the  actual  motion  better  (Der  et  al  1993). 

Slowness  measurements  across  arrays  are  notoriously  inaccurate  because  of  the  spatial 
variations  of  waveforms,  which  is  incompatible  with  the  plane  wave  propagation  model.  Figure 
1  is  a  good  indication  of  what  can  be  accomplished  with  these  two  measures  using  all  elements 
of  a  regional  array.  It  appears  that  slowness  is  still  useful  in  separating  P  and  S  type  (fast  and 
slow)  arrivals,  while  rectilinearity  is  of  little  value. 
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Figure  I.  Plot  of  recti  linearity  vs.  slowness  at  ARCESS  for  various  regional  arrivals  (plus  signs  for  Pn,  x 
for  Pg,  dot  for  Lg  and  circles  for  Sn)  for  a  set  of  events  recorded  at  ARCESS.  The  plot  shows  that  while 
rectilinearity  is  of  little  value  in  distinguishing  P  and  S  type  phases  because  of  the  large  overlap,  slowness 
separates  them  fairly  well.  Nevertheless  P  and  S  wave  types  constitute  two  mixed  groups  with  no 
separation  of  Pn  vs.  Pg  or  Sn  vs.  Lg. 


In  earlier  work  (Der  et  al  1993),  we  have  established  the  fact  that  the  polarization  patterns  of 
regional  arrivals,  unlike  simple  P  and/or  S  type  motion  patterns,  cannot  be  fully  described  by 
three  perpendicular  components  of  motion.  If  this  were  the  case,  then  it  would  be  possible  to 
derive  one  component  of  motion  from  the  other  two  with  some  simple  linear  filters,  but  such 
was  only  possible  for  Pn,  and  not  for  the  other  regional  arrivals.  Moreover,  we  have  found  that 
two,  rather  than  one,  signal  components  are  needed  to  model  Pn.  In  general,  in  array  processing 
one  must  have  fewer  independent  signal  components  than  the  number  of  sensors  in  an  array  to 
make  such  waveform  extrapolations  possible.  Sn  and  Lg  seem  to  contain  more  signal 
components  than  two. 


The  particle  motions  of  even  supposedly  simple  ‘P  type’  arrivals  seem  to  be  quite  complex  and 
three-dimensional.  The  situation  is  even  more  complex  for  other  regional  arrivals.  Some  of  the 
complexities  in  the  various  regional  arrivals  can  be  explained  by  the  fact  that  all  regional 
arrivals  are  complex  superpositions  of  multiple  reflections  in  the  crust  and  that  each  sensor 
location  is  also  associated  with  some  complex  azimuth-slowness  dependent  transfer  functions 
(often  called  site  effects)  that  distort  the  waveform.  The  latter  arise  from  geological  complexi- 
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ties  and  topography  near  each  sensor  leading  to  waveform  variations  across  even  very  small 
arrays. 

Observations  of  regional  waveforms  at  small  arrays  do  not  conform  to  the  commonly  used 
model  of  single  plane  waves  propagating  with  uniform  phase  velocities  either.  If  this  were  the 
case,  then  the  waveforms  of  regional  arrivals  should  be  identical  at  all  sensors  of  the  same  type 
and  they  are  clearly  not.  Even  if  we  are  willing  to  include  site-dependent  simple  site-transfer 
functions,  then  it  should  be  possible  to  extrapolate  the  regional  waveforms  at  one  sensor  from 
any  of  the  other  sensors.  We  have  demonstrated  ( Der  et  al  1993 )  that  this  is  not  possible  even 
for  groups  of  events  at  nearly  the  same  location. 

These  results  are  not  unexpected,  even  in  the  conventional  seismological  framework.  Regional 
arrivals  can  be  described  as  complex  superpositions  of  P  and  S  waves  multiply  reflected  in  the 
crustal  waveguide,  (Vogfjord  and  Langston  1990)  or  alternatively,  higher  modes  of  Love  and 
Rayleigh  type.  As  such  they  are  not  simple  plane  waves  (since  they  are  dispersive)  and  may 
have  complex  polarizations  even  for  laterally  homogeneous  layered  media.  The  three- 
dimensional  heterogeneity  of  the  Earth  complicates  the  matter  further,  since  this  makes  the  ‘site 
effects’,  i.e.  intersite  multichannel  transfer  functions,  azimuth  and  slowness  dependent. 

The  above  statements  should  not  be  interpreted  as  declarations  of  the  complete,  dogmatic 
repudiation  of  conventional  signal  models,  such  as  P  and  S  wave  models  and  plane  waves. 
Obviously,  they  must  have  some  validity,  (Christofferson  and  Roberts  1996)  since  we  are  able 
to  identify  P  motion  about  50%  of  the  time  by  polarization  analysis;  thus,  with  a  high  failure 
rate  and  F-K  methods,  based  on  the  assumption  of  simple  plane  waves,  still  give  valuable 
diagnostic  information  with  regards  to  slowness  and  azimuth.  We  merely  point  out  that  these 
models  do  not  describe  the  properties  of  the  observed  regional  data  well.  Moreover,  some  gross 
features  of  the  processing  schemes  we  developed  can  also  be  explained  in  terms  of  processing 
for  simpler  models,  but  not  in  detail.  If  the  simple  polarization  or  plane  wave  models  were 
exactly  valid  our  processing  schemes  would  simply  reduce  to  some  of  the  familiar  simple 
processing  schemes,  such  as  beamforming. 

In  this  report  the  term  “polarization”  describes  the  set  of  complex,  frequency  dependent  transfer 
function  relationships  among  the  various  sensors  of  a  seismic  array  that  may  have  any  mix  of 
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vertical  and  horizontal  seismometers.  The  “polarization”  changes  as  various  types  of  seismic 
phases  arrive  and  this  will  give  us  an  opportunity  to  discriminate  the  various  arrivals.  In  this 
sense,  even  an  array  containing  only  vertical  sensors  provides  polarization  information,  even 
though  most  of  the  phase  shifts  may  be  associated  with  propagation  delays.  The  observations 
mentioned  lead  to  a  type  of  generalized  multichannel  signal  models  for 


The  transfer  functions  F^oo)  are 

•  generalized  and  the  structure  is 

•  replicated  for  all  sensors  and  involves 
the  same  set  of  signals.  Plane  wave 

•  and  analytical  polarization  models 
are  special  cases  of  the  same  model. 


Figure  2.  Generalized  polarization  model.  Multiple  independent  signal  processes  are  passed  through 
generalized  transfer  functions  to  produce  a  sensor  input.  The  whole  structure  each  time  involving  different 
transfer  functions  but  the  same  set  of  signal  processes,  is  replicated  as  many  times  as  there  are  sensors. 


regional  arrivals  illustrated  in  Figure  2  assuming  that 


a)  All  of  the  regional  arrivals  contain  more  than  the  one  or  two  independent  (orthogonal) 
signal  components  assumed  to  be  present  in  P  and  S  type  arrivals  respectively. 

b)  The  waveforms  of  each  sensor  in  a  seismic  array  (including  3-component  combinations) 
are  to  be  regarded  as  outputs  of  such  multiple  signal  components  as  passed  through  some 
sets  of  multichannel  linear  filters  related  to  site-effects  and  various  frequency-dependent 
propagation  delays.  Such  multichannel  filters  will  be  assumed  to  remain  constant  for  the 
same  regional  arrivals  for  events  at  the  same  distance  and  azimuth  from  the  array. 


c)  Multiple  events  at  nearly  the  same  locations,  therefore,  will  be  assumed  to  be  independent 
replications  that  contain  the  same  set  of  multichannel  filters  and  thus  can  be  used  to  build  up 
statistics  used  in  the  detection  and  identification  of  various  regional  arrivals  at  that  distance 
and  azimuth. 
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Needless  to  say,  some  of  these  assumptions  must  be  verified  using  data.  While  such 
assumptions  are  not  useful  for  isolated  single  events,  repeated  events  at  nearly  the  same  location 
are  inevitable  consequences  of  mining  activity,  and  thus  multiple  events  at  nearly  the  same 
location  are  common.  Thus,  the  basic  idea  underlying  our  approach  is  to  replace  time  delay 
operators  and  assumed  particle  motion  phase  shifts  used  in  deterministic  signal  processing 
( Kwaerna  and  Ringdal  1992)  with  empirically  computed  sets  of  transfer  functions  and  a  general 
multi-signal  model.  Although  the  plausibility  of  a  multiple-signal  multichannel  model  has  been 
demonstrated  for  a  limited  sets  of  Kola  peninsula  events  ( Der  et  al  1993),  much  needs  to  be 
done  to  extend  these  finding  to  other  set  of  events.  In  that  paper  it  has  been  shown  that  Pn 
contained  two  signal  components  instead  of  one,  and  Sn  and  Lg  must  contain  more.  This  was 
shown  by  semi-qualitative  comparisons  of  multichannel  predictions  of  one  channel  from  the  rest 
for  a  mini-array  consisting  of  a  three-component  combination  and  a  small  tripartite  vertical 
array.  This  work  does  not  answer  the  following  questions  however: 

a)  How  many  effective  signal  components  are  present  in  Sn  and  Lg? 

b)  How  do  these  things  vary  with  distance  and  region  and  arrival? 

c)  What  statistical  framework  must  be  set  up  for  detection  and  discrimination  of  the  various 
regional  arrivals? 

d)  What  is  the  best  minimum  sensor  configuration  to  use? 

e)  What  is  the  most  parsimonious  model  (both  in  terms  of  signal  components  and  time  domain 
filter  coefficients)  that  is  still  effective? 


One  can  arrive  at  an  alternative  way  to  look  at  the  problem  by  considering  the  fact  that  a 
multisensor  seismic  array  samples  the  motion  at  a  small  portion  of  a  heterogeneous  and  possibly 
anisotropic  elastic  continuum  at  a  few  points.  Clearly,  since  these  points  are  elastically 
connected  together  they  cannot  move  independently  from  each  other.  Nevertheless,  there  must 
be,  intuitively,  some  ways  to  distinguish  the  generalized  particle  motion  of  the  various  arrivals 
as  long  as  the  source-to-array  paths  and  compositions  in  terms  of  independent  generalized  signal 
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processes  are  similar.  Conventional  analytical  particle  motion  and  propagation  models  of  wave 
types  are  of  limited  use  here. 

The  Kola  Peninsula  Data  Set 

Throughout  this  report  we  present  the  results  of  data  analyses  performed  on  a  data  set  of  Kola 
peninsula  quarry  blasts  as  recorded  at  ARCESS.  This  was  part  of  a  data  set  that  was  assembled 
at  ENSCO.  All  the  events,  quarry  blasts,  occurred  during  1990.  The  quarry  blasts  were  grouped 
into  event  sets  with  respect  to  their  locations  as  determined  by  the  IMS  (Bache  and Bratt  1990) 
and  also  considering  their  visible  characteristics  as  originating  in  areas  (quarries)  named  K1 , 

K2,  K3,  K4,  K5  and  K8.  There  is  no  independent  confirmation,  “ground  truth”,  of  the  validity  of 
these  groupings,  and  they  may  be  in  error.  In  Figure  3,  we  show  a  scatter  plot  of  their  locations 
retrieved  from  their  CMR  origin  files  together  with  the  location  determined  at  ARCESS.  These 
scatter  plots  show  some  outliers  and  overlaps  of  locations  of  events  that  were  supposed  to  have 
taken  place  in  the  same  quarry.  Note  that  two  events  supposedly  at  mine  K4  are  grouped  in 
location  with  the  cluster  of  K1,K2  and  K5  events.  One  of  the  K2  events  is  an  outlier  in  this 
location  plot.  The  K8  and  the  rest  of  the  K4  events  are  clearly  grouped  at  two  distinct  locations 
compared  to  the  rest.  The  original  association  of  these  events  with  mines  was  based  on 
waveform  characteristics  at  the  stations  of  the  Finnish  seismic  network.  A  reevaluation  of 
locations  may  be  necessary  including  the  ARCESS  data.  We  had  seven  events  at  Kl,  nine  at  K2, 
eleven  at  K4,  five  at  K5  and  five  at  K8. 
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Figure  3.  Scatter  plot  of  the  relative  locations  of  Kola  peninsula  quarry  blasts  analyzed  as  reported  by  the 
IMS.  The  blue  circles  --  K2,  red  plus  signs  --  K4,  red  x  letters  -K8. 

In  this  report  we  use  only  six-element  arrays,  subsets  of  ARCESS,  in  two  major  configurations, 
the  “mini-array”  used  in  Der  et  al  (1993)  which  consists  of  a  three-component  combination  at  AO 
combined  with  three  vertical  components  in  the  C  ring,  and  an  “all  vertical”  array  that  combines 
the  vertical  component  at  AO  with  two  vertical  components  in  the  B  ring  and  three  in  the  C  ring. 
The  choice  of  this  odd  combination  of  sensors  was  not  entirely  free,  but  dictated  by  the 
commonality  of  error  free  traces  among  events  and  mines.  We  tried  to  maximize  the  amount  of 
data  common  error  free  traces,  available.  In  processing  the  data  we  extrapolate  the  trace  AO  from 
the  rest  using  a  five-channel  multichannel  filter.  Increasing  the  number  of  traces  could  increase  the 
efficiency  of  the  processes,  provided  that  more  events  become  available.  Throughout  the  report  the 
events  are  designated  by  letter-number  combinations  such  as  K2110,  where  the  first  two 
characters  designate  the  mine  K2,  and  the  last  three  numbers  the  date  of  the  event  (in  year  1990). 
Most  of  the  events  have  high  S/N  ratios  on  the  broad-band  displays;  thus,  noise  is  not  a  major 
factor  in  our  analyses.  Table  I  contains  the  location  information  and  local  magnitudes  of  the 
events  used  for  various  data  analyses  in  this  report.  As  can  be  seen  from  the  table  the  K4  mine 
magnitudes,  as  a  whole,  are  lower  than  the  magnitudes  at  the  other  mines.  Therefore,  the  S/N 
ratios  of  these  events  are  lower  and  thus  are  less  suitable  for  designing  filters. 
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TABLE  I 


Event  ML  Lat  Long 


K1171  2.19  67.6315  33.8952 
K.1189  2.09  67.6110  33.7267 
K1224  2.18  67.5642  34.0417 
K1280  2.23  67.6881  33.8215 
K1287  2.40  67.5388  34.1983 
K.  1 301  2.43  67.6908  33.8736 
K.2054  2.92  67.6448  34.3462 
K2066  3.20  67.6148  34.0296 
K2110  2.51  67.6098  33.7326 
K.2147  2.26  67.6951  34.0329 
K2182  2.61  67.6695  34.0270 
K2219  2.06  67.3717  33.2637 
K.2246  2.30  67.5401  33.5353 
K2282  2.18  67.7470  33.7759 
K2285  2.41  67.5525  33.9008 
K4025  2.21  68.1233  32.8311 
K4118  2.03  67.9800  33.7571 
K4130  2.16  68.0203  32.9889 
K4139  2.15  68.0327  32.7990 
K4146  2.10  68.1001  33.7935 
K4178  2.06  67.5819  33.2103 
K4221  2.08  68.1549  32.4478 
K4230  2.07  68.0483  33.4415 
K4244  2.11  68.0751  33.9987 
K4270  2.10  67.6490  33.4750 
K5040  2.66  67.6357  33.7198 
K.5089  2.59  67.6388  34.1093 
K5125  3.13  67.5261  33.9092 
K5222  2.24  67.6576  34.4123 
K5278  2.70  67.6390  34.4300 
K8178  2.10  67.7142  30.9485 
K8223  2.03  67.6126  30.6305 
K8279  2.12  67.5405  30.9881 
K8286  2.17  67.5346  30.7000 
K8301  2.17  67.5552  30.6578 


8 


The  Anatomy  of  Regional  Arrivals 


In  order  to  optimize  the  processing  of  the  data  it  is  prudent  to  gain  some  understanding  of  the 
overall  properties  of  the  various  arrivals  from  the  various  quarries  as  seen  through  some  rather 
conventional  analyses.  A  simple  approach  is  to  run  data  from  some  sensors  through  a  bank  of 
band-pass  filters  to  assess  frequency  contents  and  relative  amplitudes  of  each  arrival.  In  this 
report  we  use  the  standard  Pn-Pg-Sn-Lg-Rg  designations  with  the  full  understanding  that  such 
may  not  be  very  meaningful  in  the  physical  sense.  More  careful  analyses  of  regional 
seismograms  made  clear  that  Pg  and  Lg,  for  instance,  are  likely  to  be  made  up  of  complex 
superpositions  of  various  wave  types  such  as  multiple  reflections  in  the  crust  (Vogfjord  and 
Langston  1990). 

We  have  run  some  representative  events  through  band-pass  filters.  Most  events  from  the  Kl 
and  K2  mines  (  including  the  K2  mine  events  analyzed  in  Der  et  al  1993  )  show  indications  of 
arrivals  that  can  be  tagged  as  Pn,  Pg,  Sn  and  Lg  (Figures  4  &  5).  The  Pn,  Pg  and  Sn  arrivals  are 
usually  the  most  pronounced  in  the  highest  frequency  bands,  but  Lg  is  the  strongest  in  the  3-6 
Hz  range.  At  higher  frequencies  Lg  is  mixed  up  in  the  strong  coda  of  the  three  earlier  arriving 
phases.  This  makes  the  use  of  filters  delimiting  the  Lg  in  the  3-6  Hz  band  advisable  during 
processing  of  Lg. 

The  K4  mine  events  show  relatively  weak  expressions  of  the  phases  Pg  and  Sn  at  ARCESS,  but 
have  a  very  strong  Pn  with  a  prolonged  coda  in  the  high  frequency  bands  (Figure  6).  The  Pn 
arrival  is  often,  but  not  always,  emergent.  The  Lg  phase  is  dominant  in  the  3-6  Hz  band  and,  for 
some  events,  the  Rg  in  the  lower  frequency  bands.  The  K5  mine  waveforms  are  again  similar  to 
those  from  Kl  and  K2  (Figure  7)  with  the  Pn,  Pg,  Sn  and  Lg  visible  and  distinct.  The  K8  mine 
events  have  two  diffuse  wavegroups,  P  and  S  types  with  no  individual  arrivals  seen  in  either  of 
them  (Figure  8). 

The  waveform  characteristics  show  clear  correlation  with  the  areal  distribution  of  the  events. 
Since  the  Kl,  K2  and  K5  events  are  closely  located,  their  waveform  characteristics  and 
frequency  structures  are  similar.  The  K4  events  are  apparently  distributed  over  a  wide  area,  and 
may  not  be  at  the  same  quarry  at  all.  The  K8  events  are  much  closer  to  ARCESS.  The  back- 
azimuths  of  the  K8  and  K4  events  are  also  different  from  that  of  the  K1-K2-K5  group. 
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Figure  4.  Band-pass  filtered  AO  seismograms  for  selected  K1  mine  events.  The  traces  shown  are  102 
seconds  long. 
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Figure  5.  Band-pass  filtered  AO  seismograms  for  selected  K2  mine  events.  The  traces  shown  are  102 
seconds  long. 
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Figure  6.  Band-pass  filtered  AO  seismograms  for  selected  K4  mine  events.  The  traces  shown  are  102 
seconds  long. 
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Figure  7.  Band-pass  filtered  AO  seismograms  for  selected  K5  mine  events.  The  traces  shown  are  102 
seconds  long. 
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Figure  8.  Band-pass  filtered  AO  seismograms  for  selected  K8  mine  events.  The  traces  shown  are  102 
seconds  long. 


A  common  problem  is  the  dominance  of  low  frequencies  in  the  noise,  a  common  property  of 
microseisms.  Including  these  low  frequencies  in  the  processing  initially  caused  high  correlations 
for  all  filters  designed  to  detect  various  phases  in  the  noise  windows  preceding  the  events.  When 
we  limited  the  processing  to  frequencies  above  1  Hz,  where  most  of  the  signals  dominate,  this 
problem  was  remedied. 
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OUTLINE  OF  A  THEORETICAL  FRAMEWORK  FOR  REGIONAL  ARRIVAL 
IDENTIFICATION 


General  Approach 

Assuming  that  the  multichannel  structure  of  signals  is  as  depicted  in  Figure  2,  this  would  have  the 
consequence  that  waveforms  at  any  sensor  could  be  reconstructed  from  the  waveforms  at  the  rest 
of  the  sensors  as  long  as  the  number  of  independent  signal  processes  is  less  than  the  number  of  the 
sensors  of  the  array.  The  mode  of  reconstruction  depends  on  the  matrix  of  frequency-dependent 
transfer  functions  and  not  on  detailed  time  histories  of  the  individual  processes.  As  long  as  the 
matrix  of  transfer  functions  is  the  same  for  a  set  of  events,  the  reconstruction  could  be  effectively 
accomplished  using  the  same  set  of  multichannel  filters  (Figure  9).  When  this  matrix  is  different 
few  an  event  (possibly  because  of  a  different  modal  composition,  azimuth  and  slowness  structure) 
a  noticeable  degradation  should  occur.  The  filter  design  and  performance  evaluation  can  be 
performed  in  both  the  time  and  frequency  domains  using  cross-correlation  function  or  spectral 
matrices  constructed  from  a  set  of  training  events  (quarry  blasts)  known  to  be  at  the  same 
location.  The  formulas  used  for  filter  design  are  standard  ( Shumway  1989,  Wiggins  and  Robinson 
1965).  The  time  domain  performance  evaluation  uses  the  time  domain  RMS  error  and  a  likelihood 
ratio  methodology,  and  the  frequency  domain  approach  uses  the  power  residuals  and  F  statistics. 


P  (extrapolation 
^  error) 


Filters  were  found  to  be  specific  to 
gross  path  and  arrival  type, 
but  not  specific  to  event  within  a 
group. 


Reference 

sensor 


Figure  9.  Processing  flow  of  the  small  array  polarization  processor. 
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The  time  domain  errors  or  power  residuals  should  be  small  for  events  belonging  to  the  same 
group  as  the  training  events,  and  large  for  the  events  not  belonging  to  the  group. 

Frequency  Domain  Approach 

Typically,  there  are  a  number  of  events  at  a  given  quarry  which  must  be  used  to  construct  a 
spectral  matrix  of  all  the  sensors  by  event  compounding  the  cross  products  of  the  Fourier 
transforms  of  the  time  windows  placed  over  the  various  arrivals  seen  in  the  seismogram. 

It  is  desirable  that  the  time  windows  are  placed  consistently,  although  small  inaccuracies  are  of 
no  consequence.  Typically,  we  have  lined  up  the  Pn  arrivals  at  the  sensor  AO  and  placed  the  256 
point  (6.4  second)  windows  preceding  it  by  about  50  points.  This  way  the  first  arrival  was  not 
tapered,  but  the  window  included  some  noise  that  was  rather  small  for  all  events.  The  first  and 
last  30  points  of  the  windows  were  tapered  off  using  a  cosine  taper.  The  windows  for  the  later 
phases  were  placed  at  consistent  time  intervals  relative  to  the  Pn  windows.  The  starting  times  of 
the  Pg  and  Sn  windows  again  preceded  the  visible  arrival  to  avoid  tapering  the  initial  part  of 
these  arrivals.  The  Lg  windows  were  placed  on  the  initial  large  amplitude  part  of  this  phase 
since  this  part  was  shown  to  be  more  coherent  across  arrays  ( Der  et  al  1986).  Propagation 
delays  across  the  arrays  were  not  considered  in  the  placing  the  windows,  i.e.  the  windows  for  all 
the  sensors  are  time-aligned,  since  time  delays  are  considered  implicitly  in  the  computed 
transfer  functions  and  are  small  relative  to  the  window  lengths. 


The  approach  used  in  {Der  et  al  (1993)  is  based  on  an  implicit  multichannel  model  of  the 
signals  received,  which  assumes  that  the  signals  can  be  described  as  in  terms  of  m  superposed 
independent  processes  seen  on  n  channels.  Clearly,  for  making  the  problem  tractable  it  is 
required  that  n>m,  i.e.  that  the  number  of  channels  be  larger  than  the  number  of  independent 
generalized  processes.  In  such  cases  the  model  can  be  written  in  a  partitioned  form  in  the 
frequency  domain  for  a  given  frequency 


(1) 


where  y  is  an  n  by  1  column  complex  vector  of  Fourier  components  at  frequency/,  x  is  an  m  by 
1  complex  vector  of  the  independent  signal  processes  (assumed  to  have  an  RMS  value  of  unity) 
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which  vary  from  event  to  event,  A  is  an  n  by  m  complex  matrix  if  frequency  domain  transfer 
functions,  which  would  also  incorporate  the  relative  amplitudes  of  x,  that  also  are  assumed  to  be 
nearly  the  same  for  collocated  events  seen  at  an  array.  Partitioning  the  problem  such  that  Aj  is 
an  m  by  m  matrix  and  yj  is  on  m  by  1  vector,  then  we  could,  in  principle,  design  a  multichannel 
filter  that  could  extrapolate  y2  from  the  yj  if  we  knew  A. 

Such  a  filter  would  be  obviously  A2Aj*1,  and  would  exist  as  long  as  Aj  is  not  singular.  It 
happens  practically  never  because  there  exist  some  noise  in  the  data  always.  We  could  design 
such  a  filter  if  we  knew  A,  but  obviously  we  do  not  know  either  A  or  x.  One  could  also  design  a 
filter  that  would  extrapolate  the  waveform  on  one  sensor  from  the  rest.  Such  a  filter  would  work 
for  all  events  with  the  same  A  regardless  of  the  differences  in  x. 

Nevertheless,  it  is  possible  to  design  such  multichannel  filters  from  the  data  by  using  the 
standard  frequency  domain  formulation  of  the  Wiener-Hopf  equation 

h=S~'c  (2) 

where  h  is  the  complex  multichannel  transfer  function  (an  n-1  by  1  vector)  and  S  is  the  (n- 
l)X(n-l)  complex  spectral  matrix  of  the  time  series  recorded  at  the  various  sensors,  excluding 
the  one  to  be  extrapolated,  and  c  is  the  n-1  by  1  complex  column  vector  of  the  cross  spectra 
between  the  n-1  input  processes  and  the  channel  to  be  extrapolated.  Both  S  and  c  are  obviously 
parts  of  the  total  spectral  matrix  A  of  the  array  data  and  can  be  computed  directly  from  the  array 
recordings.  As  long  the  makeup  of  the  signals  are  similar,  one  can  compute  S  by  averaging  the 
cross  products  qf  the  Fourier  transform  components  over  suites  of  n  events 

s,  =  £/*(<»)/,»  (3) 

n 

where  the  *  denotes  complex  conjugation.  In  order  to  weigh  the  various  events  with  varying 
event  magnitudes  equally,  we  normalized  the  average  absolute  amplitudes  of  each  event  to  the 
same  value  but  kept  the  relative  amplitudes  and  phases  among  sensors  the  same  for  the  same 
suite  of  events.  As  long  as  A  is  nearly  the  same  for  all  events,  the  processing  using  the  same  h 
will  perform  well  for  all  events.  In  terms  of  physics,  A  can  be  thought  of  as  a  suite  of  some 


17 


transfer  functions  that  could  be  interpreted  as  describing  the  site  responses  and  relative 
amplitudes  of  generalized  multipaths  or  multiple  signals  seen  across  the  array.  Coordinate 
rotations  and  intersensor  propagation  delays  have  been  conveniently  implicitly  absorbed  in  A 
and  needed  not  be  considered  separately.  We  are  not  in  a  position  to  extract  and  interpret  either 
A  or  x  in  terms  of  transfer  functions  or  physical  signal  processes,  but  may  still  exploit  the 
eventual  similarities  in  the  structures  of  quarry  blast  signals  arriving  along  near  identical  paths. 
There  is  no  guarantee  that  A  remains  nearly  constant  for  such  suites  of  collocated  events,  but 
apparently  this  seems  to  be  the  case  in  most  cases  as  the  results  below  and  in  the  original  paper 
{Der  et  al  1993)  show.  The  method  rests  on  the  homogeneity  of  transfer  functions  that  has  to  be 
tested  (Appendix  A). 

Typically,  it  is  advantageous  to  stabilize  the  inversion  of  S  by  regularizing  it,  i.e.  adding  to  the 
diagonal  a  positive  multiple  of  the  unit  matrix.  Two  kinds  of  regularizing  was  implemented  in 
the  code.  The  first  consists  of  multiplying  the  diagonal  elements  of  the  spectral  matrix  A  with 
some  factor,  the  second  is  by  adding  a  true  unit  matrix  which  was  a  set  multiple  of  the 
maximum  of  the  diagonal  of  A.  We  have  found  no  visible  differences  in  the  results  of  the  two. 
We  have  applied  the  first  most  of  the  time  with  a  multiplying  factor  set  at  1 .5. 

In  the  case  where  the  seismic  noise  is  a  problem,  in  doing  the  waveform  extrapolation,  the 
Wiener  filters  can  be  modified  by  including  the  noise  spectral  matrix  N  in  the  design 

h  =  (s+n)”'  c.  (4) 

In  order  to  do  this  one  needs  to  estimate  N  from  the  seismogram  preceding  each  individual 
event  to  be  processed.  This  way  the  process  can  be  optimized  for  the  given  background  noise 
field.  Obviously  it  is  advantageous  from  the  reliability  point  of  view  to  estimate  S  and  c 
separately  from  high  S/N  events  and  design  h  from  the  formula  prior  to  applying  it  to  any  low 
S/N  events. 

Time  Domain  Approach 

The  scheme  of  identifying  a  regional  arrival  type  from  the  comparison  of  an  extrapolated 
sensor  output  to  an  actual  one  can  be  reduced  to  a  statistical  testing  of  a  set  of  residuals  to 


18 


determine  which  one  of  the  regional  phase  models,  defined  in  terms  of  multichannel  filters,  is 
the  most  likely  one.  This  is  similar  to  the  schemes  described  by  Basseville  and  Nikiforov  (1993) 
to  determine  which  of  a  set  ARMA  models  fits  a  set  of  data  best.  The  residuals  of  extrapolation 
of  one  output  from  the  rest  can  be  written  as 

l=N 

£k,i  =  Sk,i  IX,-,  / '/./  (^) 

jtk  l=*-N 

where  fjj  are  the  filter  weights  applied  to  the  the  traces  Sjj  with  the  first  sum  over  all  sensors 
excluding  the  sensor  output  k  to  be  extrapolated  from  the  rest.  The  length  of  the  time-symmetric 
filter  is  2N+1.  The  weights  can  be  obtained  by  solving  the  well-known  time  domain 
multichannel  filter  design  equations  ( Wiggins  and  Robinson  1963),  where  the  r/  are 
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block  kxk  autocorrelation  matrices,  subscripts  are  for  each  time  lag  /,  where  k  is  the  number  of 
input  traces,  and  the  gj  are  the  cross  correlation  vectors  between  the  inputs  and  the  trace 
designated  as  the  desired  output,  and  fi  are  filter  weights.  The  large  matrix  is  block-Toeplitz  in 
form  and  is  thus  amenable  to  iterative  solution  methods.  These  can  be  obtained  by  ensemble¬ 
averaging  the  corresponding  correlation  functions  over  suites  of  collocated  events  for  each 
specific  arrival.  The  weights  will  be  specific  to  each  arrival  the  sample  windows  will  be 
centered  on.  We  have  implemented  the  recursive  algorithm  of  Wiggins  and  Robinson  (1965)  for 
use  in  our  work. 

The  correlation  functions  were  computed  in  a  manner  similar  to  the  spectral  matrices  above  by 
summing  the  cross-correlation  functions  among  sensor  pairs  of  the  various  events  after  RMS 
event-normalization  of  the  waveforms  to  avoid  the  domination  of  the  cross-correlation  function 
matrices  by  the  events  with  the  largest  magnitudes. 
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DATA  ANALYSIS 


Frequency  Domain. 

The  time  domain  method  as  described  above  was  applied  to  various  data  sets.  The  “mini-array’ 
consists  of  a  three-component  set  at  AO  site  of  ARCESS  combined  with  three  vertical  sensors  in 
the  C  ring  (Figure  10).  Applying  the  multichannel  waveform  equalization  method  to  the  two  sets 
of  data  results  in  good  replicas  of  the  waveforms  of  all  phases.  Figures  11  to  14  show  the  results 
for  some  examples  of  Pn,  Pg,  Sn  and  Lg  for  the  original  K2  data  set. 


ICiOm 


Figure  10.  ARCESS  array  configuration. 
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Figure  1 1 .  Examples  of  Pn  waveform  extrapolations  using  the  “mini-array”  and  the  leave-one-out  method. 
All  of  the  waveform  matches  between  the  original  AO  waveforms  (top)  and  the  waveform  extrapolations 
(below)  are  excellent.  Event  identifications  and  correlation  coefficients  are  plotted  above  each  pair  of 
waveforms.  The  waveform  sections  shown  are  6.4  seconds  long. 
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Figure  12.  Pg  waveform  extrapolations  using  the  “mini-array”  and  the  leave-one-out  method.  All  of  the 
waveform  matches  between  the  original  AO  waveforms  (top)  and  the  waveform  extrapolations  (below)  are 
good.  Event  identifications  and  correlation  coefficients  are  plotted  above  each  pair  of  waveforms.  The 
waveform  sections  shown  are  6.4  seconds  long. 


22 


si : kXliS4.dat 

V  H'-JhMhfr-Ut 

.  .  . 


A  :  W. !  tfJ.Hftf 

ref  •  V.IWTSKE-OI 

i:*.d . . 


u.f  6  .  01 

raid . ,  , 


Figure  13.  Sn  waveform  extrapolations  using  the  “mini-array”  and  the  leave-one-out  method.  The 
waveform  matches  between  the  original  AO  waveforms  (top)  and  the  waveform  extrapolations  (below)  are 
reasonably  good.  Event  identifications  and  correlation  coefficients  are  plotted  above  each  pair  of 
waveforms.  The  waveform  sections  shown  are  6.4  seconds  long. 
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Figure  14.  Lg  waveform  extrapolations  using  the  “mini-array”  and  the  leave-one-out  method.  For  this 
arrival  The  waveform  matches  between  the  original  AO  waveforms  (top)  and  the  waveform  extrapolations 
(below)  are  still  good  but  for  a  few  events  the  correlation  coefficients  are  low,  such  as  the  one  pair  at  the 
bottom.  Event  identifications  and  correlation  coefficients  are  plotted  above  each  pair  of  waveforms.  The 
waveform  sections  shown  are  6.4  seconds  long. 
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The  complete  runs  on  events  were  performed  the  following  way.  The  frequency  domain  filters 
were  computed  from  the  event-compounded  spectral  matrices  of  a  learning  data  set.  These 
filters  were  applied,  in  the  frequency  domain,  over  half-overlapping  256  point  windows  over  the 
whole  recording  to  construct  the  extrapolated  reference  trace  arbitrarily  chosen  to  be  the 
vertical  at  site  AO.  Normalized  correlation  coefficients  were  then  plotted  as  a  function  of  time 
for  the  four  filter  outputs  designed  to  extrapolate  the  wavefields  appropriate  to  Pn,  Pg,  Sn  and 
Lg.  Note  that  the  information  makes  no  use  at  all  of  the  variations  in  the  absolute  amplitudes  of 
the  input  which  define  a  regional  arrival  to  an  analyst.  Neither  does  the  increase  of  correlations 
at  the  arrival  of  Pn  has  anything  to  do  with  the  increase  of  wave  amplitudes  at  that  time.  This 
increase  of  correlation  coefficients  would  occur  even  if  the  amplitudes  did  not  change.  It  merely 
shows  that  the  motion  is  more  “Pn-like”  in  the  windows  than  in  the  windows  preceding  it. 

Inspecting  the  results  of  complete  processing  runs  on  some  complete  events  show  that,  in 
general,  it  is  possible  to  distinguish  between  Pn  and  Pg,  and  Sn  and  Lg  in  most  cases.  There  is 
no  ambiguity  in  telling  the  P  type  and  S  type  phases  apart.  While  Pn  is  the  first  maximum  on  the 
Pn  runs,  Pg  may  be  still  smaller  that  Pn  in  the  Pg  runs,  but  it  is  generally  more  pronounced  and 
enhanced  relative  to  Pn  in  the  Pg  run.  Figure  15  shows  results  for  some  events  in  the  K2  mine. 
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Event  K2285 


Figure  15.  Processed  events  from  the  K2  mine  using  the  “mini-array”  configuration.  The  seismogram  of 
the  AO  sensor  is  shown  on  top.  The  correlation  coefficient  between  256  point  (6.4  sec)  half-overlapping 
windows  on  the  AO  trace  and  the  extrapolated  waveform  is  plotted  as  a  function  of  time  along  the  time 
axis.  The  vertical  line  segments  on  the  left  indicate  the  0-1.0  range  in  the  size  of  correlation  coefficients. 
The  results  for  filter  sets  designed  to  identify  Pn,  Pg,  Sn  and  Lg  are  plotted  from  top  to  bottom.  Note  that 
the  maxima  tend  to  occur  at  the  right  arrivals  (v  marks).  Often  the  Pg  (or  other  phase)  is  identified  only  by 
a  large  increase,  not  an  absolute  maximum,  at  the  time  of  the  arrivals  o  (?  Marks).  Nevertheless,  such 
increases  often  identify  the  phase  well.  The  total  length  of  the  records  processed  in  this  type  of  figure  was 
102.4  seconds  throughout  this  report. 
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Event  K.2285 


Figure  16.  K2  events  processed  by  the  leave-one-out  technique.  In  these  runs  we  have  excluded 
frequencies  below  1  Hz,  where  the  background  noise  dominates  and/or  it  is  similar  to  the  Pn  because  of 
the  longer  wavelengths.  This  exclusion  from  the  processing  was  effected  by  not  using  these  low 
frequencies  in  the  waveform  reconstruction.  The  resulting  small  loss  in  the  size  of  correlation  coefficient 
apparently  was  not  serious.  Length  of  these  time  traces  is  1 02.4  seconds. 


An  important  question,  since  we  are  doing  no  spectral  smoothing  is,  how  many  events  are 
needed  for  the  scheme  to  work.  In  order  to  test  the  robustness  of  the  scheme,  we  have  made  runs 
with  the  K.2  data  set  where  the  event  being  processed  was  excluded  from  the  training  data  set 
used  for  computing  the  filters.  These  processing  runs  still  identified  the  appropriate  arrivals 
quite  well  (Figure  16)  with  some  exceptions  for  Lg.  In  general,  Pn  waveform  extrapolations 
were  quite  robust  when  the  leave-one-out  method  was  applied,  but  the  robustness  decreased 
from  Pn  to  Lg.  It  appears  that  one  should  not  use  much  less  than  nine  events  as  a  learning  set. 
These  should  be  in  the  same  general  area,  but  need  not  be  in  the  same  quarry  ( Der  et  al  1993). 

It  is  interesting  to  inspect  the  time  domain  equivalents,  impulse  responses,  of  the  filters  being 
used  for  identifying  Pn,  Pg,  Sn  and  Lg.  In  Figures  17-20  we  show  these  for  the  K2  mine.  These 
show  that  in  a  sense  what  we  are  doing  is  similar  to  beamforming.  The  highest  amplitude 
portions  of  the  impulse  responses  follow  approximately  the  propagation  time  delays 
appropriate  to  the  phases.  Note  that  these  delays  increase  from  Pn  to  Lg.  Nevertheless,  in  pure 
beamforming  we  should  have  only  delayed  delta  functions  with  zeroes  in  between,  while  our 
responses  are  quite  complex,  and  thus  we  are  doing  a  complex  filter-sum  instead  of  delay  sum. 
It  is  interesting  to  note  that  for  Pn  reconstruction,  we  now  use  mostly  the  available  C  ring 
vertical  sensors,  largely  ignoring  the  horizontal  components,  while  in  Paper  I  we  have  shown 
that  the  two  horizontal  components  of  motion  at  AO  can  accomplish  the  same  goal  quite  well. 
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When  regularizing  the  spectral  matrices  we  have  constrained  the  processor  such  that  it  does  the 
least  amount  of  waveform  transformation  to  accomplish  the  goal  of  effective  wavefield 
extrapolation. 
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Figure  17.  Impulse  responses  of  filters  designed  to  extrapolate  the  Pn  wavefield  for  events  at  the  K2  mine 
using  the  “mini-array”  from  ARCESS.  The  zero-lag  time  is  at  the  middle  of  each  waveform.  Time  length 
of  these  impulse  responses  is  6.4  seconds. 
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Figure  18.  Impulse  responses  of  filters  designed  to  extrapolate  the  Pg  wavefield  for  events  at  the  K2  mine 
using  the  “mini-array”  from  ARCESS.  Time  length  of  these  impulse  responses  is  6.4  seconds. 


30 


AOsn 

AOse 


Figure  19.  Impulse  responses  of  filters  designed  to  extrapolate  the  Sn  wavefield  for  events  at  the  K2  mine 
using  the  “mini-array”  from  ARCESS.  Time  length  of  these  impulse  responses  is  6.4  seconds. 
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Figure  20.  Impulse  responses  of  filters  designed  to  extrapolate  the  Lg  wavefield  for  events  at  the  K2  mine 
using  the  “mini-array”  from  ARCESS.  Time  length  of  these  impulse  responses  is  6.4  seconds. 


As  an  alternative  sensor  configuration  we  have  also  evaluated  a  combination  of  three  sensors  of 
the  C  ring,  the  vertical  sensor  at  AO  and  two  vertical  sensors  in  the  B  ring,  (many  of  these 
choices  were  dictated  by  the  availability  of  data  channels  without  errors).  We  have  recovered 
data  for  this  combination  for  all  the  Kola  mines  except  K3.  This  configuration  has  the  drawback 
of  small  effective  aperture,  but  it  discards  the  two  horizontal  components  at  the  center  that  are 
not  effectively  utilized  (as  revealed  by  Figures  17-20)  for  the  later  arrivals  in  the  “mini-array”, 
especially  Sn  and  Lg.  In  retrospect,  a  better  configuration  could  consist  of  the  vertical  sensors 
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of  the  whole  C  ring  combined  with  the  vertical  sensor  at  AO. 


The  runs  for  this  configuration  again  show  the  effective  identification  of  the  four  major  arrivals 
(Figure  21)  for  the  K2  mine,  but  the  differences  among  them  are  less  pronounced.  The 
beamforming  has  less  discrimination  power,  because  of  the  closeness  of  AO  to  the  B  ring 
sensors.  Nevertheless,  in  such  processing  although  we  are  giving  more  weight  to  the  closest 
sensors  in  extrapolating  a  sensor  output  from  the  rest,  this  is  not  simple  copying-over  with  a 
time  delay.  Actually,  it  can  be  verified  by  inspection  of  the  extrapolated  traces  that  they  are 
better  reproductions  of  the  AO  waveform  than  the  waveforms  seen  at  the  nearest  B  ring  sensors. 


Event  K.2054 


Event K22 19 
Figure  2 1 . 
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Event  K2285 

Figure  21.  Processed  events  from  the  K2  mine  using  the  all-verticals  configuration. 


The  impulse  response  of  the  filters  show  that  the  B  ring  sensors,  which  have  waveforms  closest 
to  that  of  AO,  are  used  with  the  greatest  weights  (Figures  22  to  25).  The  relative  weighting  and 
the  filter  waveforms  are  quite  complex 
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Figure  22.  Impulse  responses  of  filters  designed  to  extrapolate  the  Pn  wavefield  for  events  at  the  K2  mine 
using  the  all-vertical  array  from  ARCESS.  The  zero-lag  time  is  in  the  middle  of  each  waveform.  Time 
length  of  these  impulse  responses  is  6.4  seconds. 
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Figure  23.  Impulse  responses  of  filters  designed  to  extrapolate  the  Pg  wavefield  for  events  at  the  K2  mine 
using  the  all  -vertical  array  from  ARCESS.  The  zero-lag  time  is  in  the  middle  of  each  waveform.  Time 
length  of  these  impulse  responses  is  6.4  seconds. 
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Figure  24.  Impulse  responses  of  filters  designed  to  extrapolate  the  Sn  wavefield  for  events  at  the  K2  mine 
using  the  all-vertical  array  from  ARCESS.  The  zero-lag  time  is  in  the  middle  of  each  waveform.  Time 
length  of  these  impulse  responses  is  6.4  seconds. 
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Figure  25.  Impulse  responses  of  filters  designed  to  extrapolate  the  Lg  wavefield  for  events  at  the  K2  mine 
using  the  all-vertical  array  from  ARCESS.  The  zero-lag  time  is  in  the  middle  of  each  waveform.  Time 
length  of  these  impulse  responses  is  6.4  seconds. 


We  have  done  some  tests  in  using  the  filters  designed  for  one  mine  (using  the  all-vertical  array 
in  all  these  tests)  on  events  of  the  other  mines.  Remembering  the  reported  locations  of  the 
events  the  K.1,  K.2  and  K5  mines  appear  to  be  located  closer  than  the  K.4  and  K8  mines  (if  these 
locations  can  be  believed).  When  the  filter  set  of  K.2  is  applied  to  the  events  of  K4  mine  the 
correlations  are  lower  and  it  is  no  longer  possible  to  see  any  differences  in  the  Pn  and  Lg  pair  or 
the  Sn  and  Lg  pair  (Figure  26).  Nevertheless,  P  and  S  type  phases  still  produce  maxima  in 
different  parts  of  the  seismogram.  This  was  expected  based  on  the  different  azimuth  and 
distance  of  the  K4  mine. 
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On  the  other  hand,  applying  the  filters  designed  for  the  K2  mine  to  events  in  the  K.1  and  K5 
mines  the  phase  identification  was  not  entirely  unsuccessful  (Figures  27  and  28),  which  was 
also  expected  based  on  their  similar  locations. 
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Figure  27.  K.1  mine  events  processed  with  the  filter  set  designed  from  K2  mine  data. 
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Figure  29.  K2  filters  applied  to  K8  mine  events. 


On  the  other  hand,  application  of  the  K2  filters  to  the  events  at  K8  is  again  a  total  failure 
(Figure  29).  We  still  are  able  to  tell  the  P  and  S  type  phases  apart,  but  there  is  no  Pg  or  Sn,  and 
they  are  also  unclear  in  the  original  seismograms. 
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Time  Domain 


Much  of  the  work  was  spent  on  developing  a  frequency-domain  implementation  of  the  regional 
phase  identification  method  of  Der  et  al  (1993).  Thus,  we  transformed  the  frequency  domain 
filters  that  equalize  the  waveforms  into  the  time  domain.  These  will  be  applied  as  convolution 
filters  in  a  continuous  fashion.  A  drawback  of  this  approach  is  that  these  filters  are  non-causal 
and  introduce  precursors.  Nevertheless,  since  our  purpose  is  phase  identification  and  the 
precursors  have  very  little  energy,  this  is  not  a  problem. 

As  described  by  Der  et  al  (1993),  we  design  digital  filters  for  predicting,  or  more  correctly, 
extrapolating  the  center  vertical  seismometer  waveform  from  the  rest  of  the  mini-array  (two 
collocated  horizontals  and  a  tripartite  configuration  of  vertical  from  the  C  ring)  based  on  the 
event  compounded  spectral  matrix  derived  from  a  set  of  events.  The  filters  then  can  be  applied 
to  individual,  collocated  events,  even  those  not  including  the  learning  set.  The  success  of  the 
waveform  extrapolation  identifies  the  phase  on  which  the  filter  was  designed.  We  have 
concentrated  on  the  Pn  and  Lg  phases  which  we  regard  as  the  key  phases  in  locating  an  event. 
The  same  approach  can,  and  has  been,  applied  to  other  identifiable  regional  arrivals  as  well, 
such  as  Sn  and  Pg  (Der  et  al  1993). 

We  show  the  time  domain  filters  for  the  extrapolation  of  Pn  phase  waveforms  in  Figure  30A. 
These  were  designed  by  inverse  Fourier  transforming  the  frequency  domain  filters  in  Der  et  al 
(1993)  and  tapering  off  the  ends.  These  filters  are  quite  complicated,  but  the  relative  amplitudes 
show  that  it  is  mostly  the  two  horizontal  components  that  contribute  to  the  prediction.  The 
forms  of  filters  do  not  lend  themselves  to  easy  intuitive  interpretation  since  the  observed 
waveforms,  being  influenced  by  distortion  due  to  local  geology,  are  quite  complex  themselves. 
Nevertheless,  the  details  of  original  and  predicted  Pn  waveforms  (Figure  30B  )  are  quite 
similar,  and  their  cross-products,  which  would  contribute  to  the  cross-correlation  coefficients, 
are  predominantly  positive  in  polarity  (lowermost  traces).  Figure  30C  shows  the  original  AO 
vertical  compared  to  the  prediction  and  the  smoothed  cross-product  for  the  whole  recording  of 
the  same  events.  In  this  display  no  normalization  with  respect  to  the  power  in  the  traces  of  the 
cross  product  has  been  performed  and  the  relative  amplitudes  of  various  arrivals  still  play  a  role. 
Nevertheless,  it  is  clear  that  the  processing  enhances  Pn  relative  to  the  other  phases,  even  in  the 
cases  (not  shown)  where  Lg  was  much  stronger  than  Pn. 
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Similarly,  the  filters  designed  to  process  the  Lg  phases  are  also  much  more  complex  than  one 
would  expect  for  simple  beamforming,  since  such  filters  would  contain  only  a  few  spikes  at  the 
correct  time  lags.  In  this  case  the  main  contributions  come  from  the  rest  of  vertical  traces 
(Figure  31A).  The  original  and  predicted  Lg  waveforms  are  quite  similar  again  in  many  details 
(Figure  31B)  and  their  cross-products  are  again  predominantly  positive.  The  waveform 
matching  procedure  did  not  work  as  well  for  Lg  as  for  Pn.  Figures  31C  show  the  original  AO 
vertical  compared  to  the  prediction  of  Lg  and  the  smoothed  cross-product  for  the  whole  AO 
sensor  recordings  of  the  same  events.  It  is  clear  that  the  processing  enhanced  Lg  considerably 
relative  to  the  other  phases.  Not  all  of  the  processed  events  show  enhancement  of  Lg  only  as 
shown  here,  however,  some  enhance  Sn  as  well.  This  can  be  expected,  since  both  Sn  and  Lg 
contain  considerable  SV  components. 

The  performance  of  filters  designed  entirely  in  the  time  domain  using  correlation  functions  and 
the  recursive  block  Toeplitz  matrix  inversion  procedure  (Wiggins  and  Robinson  1965)  was 
comparable  to  the  results  presented  above  for  Pn  even  when  very  short  31  point  filters  were 
used.  The  advantage  of  such  approach  is  that  the  filters  thus  designed  are  optimum  for  their  time 
length,  instead  of  being  truncated  as  above.  The  finding  of  the  most  effective  and  parsimonious 
filter  lengths  is  the  subject  of  our  ongoing  work.  For  Lg  the  short  filters  did  not  work  well,  we 
are  presently  modifying  and  testing  the  program  to  compute  longer  time  domain  filters,  but 
shorter  than  the  160  point  filters  shown  above.  The  need  for  testing  this  is  dictated  by  the  fact 
that  the  computation  of  longer  filters  may  make  the  recursive  algorithm  become  unstable. 
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Figure  30.  A)  Time  domain  filters  for  the  enhancement  of  Pn,  length  160  points=4  seconds.  B)  Detail  of 
AO  Pn  waveform  for  event  1990219  (top)  vs.  the  predicted  one  (middle)  and  their  cross  product  (bottom). 
Time  length  6.4  seconds.  C)  Processed  trace  (top),  actual  seismogram  (middle)  and  the  leaky-integrated 
cross  product.  Time  length  102.4  seconds. 
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Figure  31.  A)  Time  domain  filters  for  the  enhancement  of  Lg,  length  160  points=4  seconds.  B)  Detail  of 
AO  Lg  waveform  for  event  19901 10  (top)  vs.  the  predicted  one  (middle)  and  their  cross  product  (bottom). 
Time  length  6.4  seconds.  C)  Processed  trace  (top),  actual  seismogram  (middle)  and  the  leaky-integrated 
cross  product.  Time  length  102.4  seconds. 
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SUMMARY  OF  POLARIZATION  PROCESSING  RESULTS 


The  generalized  polarization  method  for  automatically  identifying  regional  arrivals  appears  to 
be  successful  for  groups  of  events  that  are  close  in  their  back  azimuths  and  slownesses  of  their 
phases  at  a  small  array,  but  not  necessarily  closely  spaced  in  the  absolute  sense.  Such  event 
groups  may  be  actually  scattered  over  distances  of  several  tens  of  kilometers  when  the  array  to 
event  distances  are  of  the  order  of  350  km.  Gross  features  of  the  generalized  polarization  filter 
operators  can  be  explained  in  terms  of  the  conventional  concepts  such  as  slownesses  and  body 
wave  polarizations,  but  in  more  detail  they  are  quite  different.  The  technique  breaks  down  as  the 
events  scatter  over  larger  areas.  The  methods  seems  to  remain  stable  when  the  event  to  be 
processed  is  left  out  from  the  suite  used  to  design  the  processor.  The  efficacy  of  the  process 
seems  to  degrade,  however,  from  Pn  to  Lg.  This  is  probably  caused  by  the  greater  complexity, 
in  terms  of  modes  or  signal  processes,  of  the  later  arrivals. 
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PHASE  ARRIVAL  TIME  ESTIMATION  AT  REGIONAL  DISTANCES  USING  THE 
CUSUM  ALGORITHM 

Introduction 

Automatic  phase  arrival  time  estimation  is  of  considerable  interest  because  of  the  need  for  rapid 
location  and  identification  of  numerous  seismic  events  by  networks  monitoring  natural  and 
man-made  seismic  activity.  Times  of  seismic  “phase”  arrivals  are  defined  as  times  where  some 
visible  characteristic,  such  as  amplitude,  frequency  content  or  wave  polarization  changed  in 
some  recording.  Typically,  regional  arrivals  are  high  frequency,  broadband,  emergent  wave 
groups  containing  numerous  cycles.  Later  arrivals  generally  have  no  clear,  impulsive 
waveforms  and  are  preceded  by  the  codas  of  earlier  ones,  and  their  onset  times  can  only  be 
defined  to  within  a  few  cycles. 

A  human  operator  often  can  spot  arrivals  by  noticing  subtle  changes  in  the  frequency  content  in 
noisy  records.  On  the  other  hand,  when  such  human  capabilities  are  needed,  it  simply  indicates  a 
failure  of  applying  appropriate  frequency  domain  prefiltering  that  could  have  already  produced 
an  enhanced  amplitude  contrast  between  the  noise  and  the  arrival.  Once  the  amplitude  contrast 
is  enhanced,  the  contrast  in  frequency  content  will  be  much  less  noticeable.  Appropriate 
prefiltering  in  frequency  is  a  prerequisite  for  further  onset  time  determination  regardless 
whether  it  is  manual  or  automatic.  Various  schemes  for  optimum  pre-filtering  were  suggested 
and  most  of  these  seem  to  work  well.  Kvaema’s  definition  of  “usable  bandwidth”  ( Kvaerna 
1995,  1996a)  and  S/N^  (in  amplitude)  or  noise-adaptive  predictive  filtering  enhances  the 
amplitude  contrast  between  the  noise  background  preceding  the  onset  of  the  signal  at  the 
expense  of  decreasing  the  visible  frequency  contrast  between  the  signal  and  noise.  Thus, 
according  to  experience  in  noisy  signals  thus  processed  it  is  hard  to  spot  the  signal  arrivals  by 
the  changing  frequency  content.  The  biggest  visual  differences  between  the  signal  and  noise 
parts  of  raw  regional  seismograms  are  due  to  the  presence  of  low  frequency  microseisms,  which 
are  eliminated  by  the  prefiltering  processes.  To  illustrate  this  point,  in  Figure  32  we  show 
examples  of  S/N^  filtered  signals  and  noise  and  in  Figure  33  simple  2.5-14  Hz  band-pass 
filtered  signals  which  both  emphasize  the  usable  bandwidth.  None  of  these  examples  show  any 
noticeable  changes  in  the  frequency  content  upon  the  arrival  of  the  signals,  and  the  signal 
arrivals  are  noticed  only  because  of  the  changes  in  the  average  trace  amplitudes. 
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Automated  equivalents  of  subjective  human  operations  requires  the  definition  of  some  statistics 
(e.g.  amplitude,  AR  model,  polarization  parameters)  and  some  algorithm  that  can  define  the 
time  when  such  parameters  changed.  The  standard  short-term  average  over  long-term  average 
ratio  (STA/LTA)  algorithm  is  well  suited  for  arrival  detection  but  not  for  precise  arrival  time 
estimation  because  of  the  long  delay  associated  with  any  significant  change  in  the  STA.  Most 
methods  proposed  for  precise  arrival  time  estimation  typically  attempt  to  detect  the  first  point 
where  some  statistic  suddenly  changes.  Many  of  these  may  work  for  most  first  arriving  Pn  if  the 
signal  is  impulsive,  but  do  poorly  with  emergent  phases  or  secondary  arrivals  typically  seen  at 
regional  distances.  The  situation  is  more  difficult  in  picking  the  later  arrivals  in  regional 
seismograms.  In  that  case  there  is  little  or  no  frequency  contrast.  Typically,  the  spectra  of  the 
various  regional  arrivals  are  quite  similar,  for  quarry  blast  sources  they  are  nearly  identical 
(Baumgardt  and  Ziegler  1986).  There  are,  thus,  no  spectral  differences  to  exploit  and  the 
determination  on  onset  times  must  depend  mostly  on  amplitude  or  polarization  changes." 
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There  is  a  basic  philosophical  issue  associated  with  the  definition  of  regional  “arrival”  onset 
times.  Typically,  most  current  routine  analyst-guided  processing  is  based  on  the  picture  of 
regional  seismograms  being  made  up  of  a  few  phases  such  as  Pn,  Pg  ,  Sn  and  Lg  typical  of 
Scandinavia  and  Europe.  Analysts  are  trained  to  recognize  and  pick  these  particular  arrivals. 
Numerous  studies  have  shown,  however,  that  regional  seismograms  are  much  more  complex 
and  could,  more  appropriately,  be  described  in  terms  of  various  physically  meaningful  wave 
groups  such  as  various  guided  wave  groups  and  reflections  such  as  multiple  PmP  and  SmS 
( Vogfford  and  Langston  1990).  An  automatic  processor  will  thus  find  more  arrivals  than  an 
analyst  indoctrinated  in  the  Pn-Pg-Sn-Lg  framework  and  these  may  not  exactly  correspond.  This 
will  be  even  more  true  for  seismograms  recorded  in  tectonic-geologic  environments  different 
from  Scandinavia  which  may  be  quite  different.  In  cases  of  such  conflicts  the  automatic 
processor  may  often  be  right  and  the  analyst  wrong.  In  this  report,  we  still  use  the  Pn-Pg-Sn-Lg 
framework  only  because  it  seems  to  fit  the  data,  not  because  we  have  a  faith  in  the  absolute 
global  validity  of  it. 

Accurate  onset  time  determination  is  important  because  the  accuracy  of  regional  location 
depends  on  errors  in  the  onset  time  picks.  If  the  onset  of  Pn  is  well  above  the  background  noise 
level  (after  appropriate  frequency  filtering),  the  associated  onset  times  should  be  used  in 
location.  On  the  other  hand,  Pn  often  has  lower  amplitudes  than  some  of  the  later  arrivals.  In 
such  cases  later  arrivals  must  be  utilized.  Moreover,  the  ability  to  determine  onset  times  of  later 
arrivals  automatically  is  useful  for  the  following  reasons:  first,  it  helps  to  identify  the  nature 
and  approximate  distance  of  the  event  by  the  pattern  of  envelope  amplitudes  and  onset  times, 
second,  some  of  the  regional  spectral  discriminants  based  on  spectral  ratios  of  various  wave 
groups  can  be  applied  automatically. 

A  related  issue  of  practical  importance  is  that  of  the  late  onset  time  picks  for  emergent  signals 
buried  in  high  background  noise.  Many  body  wave  arrivals  start  with  a  small  precursor  or  build 
up  in  amplitude  gradually.  In  such  cases  the  beginning  of  the  signal  is  missed  and  the  onset  time 
estimate  will  be  late  ( Kvaerna  1996a, b,  Douglas  et  al  1997).  In  our  opinion  this  situation  has 
no  remedy,  because,  as  we  shall  show  below,  both  the  analyst  and  any  automatic  processor  will 
pick  the  onset  time  late. 
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The  CUSUM  algorithms  described  in  Basseville  and  Nikiforov  (1993)  were  designed  for 
pinpointing  the  time  of  a  change  in  a  system  and  have  their  primary  applications  in  quality 
control  and  machine  diagnostics.  The  basic  idea  is  detecting  changes  in  the  trends  of  the 
cumulative  sum  of  some  suitable  statistic  that  abruptly  changes  with  time  as  the  properties  of 
the  time  series  change.  It  is  much  easier  to  see  and  quantify  a  change  in  a  trend,  than  to  pinpoint 
the  exact  time  of  the  first  point  where  the  change  occurred.  The  most  appropriate  statistic  to  be 
used  depends  on  the  particular  change  to  be  detected  and  will  be  varied.  In  the  case  of  seismic 
arrival  onset  determination  case  we  look  for  changes  in  signal  amplitudes,  spectral  contents  (AR 
model  residuals)  and  polarization  (cross-extrapolation  residuals)  or  any  combinations  of  these. 
The  CUSUM-based  methods  have  indeed  been  used  for  determining  P  onset  times  in  the  past 
( Nikiforov  and  Tikhonov  1986,  Nikiforov  et  al  1989),  but  not  in  a  combination  with  simulated 
annealing  as  below. 

A  simple  application  of  the  method  is  to  compute  the  cumulative  sum  of  the  chosen  statistic  and 
subtract  a  linear  ramp  from  it.  This  way  the  changes  in  the  trend  of  the  CUSUM  are  converted 
to  minima.  Numerous  methods  for  automatic  finding  of  minima  for  empirical  functions  exist. 
The  simplest  is  to  find  absolute  minima  in  the  time  windows  where  the  arrival  times  of  phases 
were  expected.  Such  time  windows  can  be  provided  by  the  generalized  polarization  analysis 
method  described  in  the  first  part  of  this  report.  This  method,  although  crude,  seems  to  be  quite 
successful  (Figures  35  and  36).  These  are  results  for  the  Khibiny  data  set  ( Mykkeltveit ,  S. 
1992). 


50 


NMIMIMHI 

Mga 

21299 


L  :  .  . 

PE 

r"  - - - — 

^ 

- - ^ 

itinnnM 

22242 


A0" 

1 

iA _ _ _ 

1 

PE 

□ 

csuK 

r — 

! — j 

K  " 

10  t«o ©not 

21293 


Mm 

m 

p"i»r 

csuSk 

fFfrT 

n 

to  Eaeondm 

22264 


Figure  35.  Application  of  the  CUSUM  method  to  four  events  in  the  Khibiny  data  set.  The  top  trace  is  the 
original  single  sensor  output,  below  it  is  the  adaptive  AR  filtered  version.  At  the  bottom  is  the  cumulative 
sum  of  the  absolute  values  of  the  second  trace  with  the  time  picks  (vertical  lines)  for  Pn,  Pg,  Sn  and  Lg. 
Note  that  in  spite  of  the  ill-defined  beginnings  of  each  phase  the  minima  are  quite  well  defined.  Absolute 
minima  in  certain  time  frames  were  picked.  These  examples  of  onset  time  estimation  are  based  on 
amplitude  changes  of  course. 


Figure  36.  Application  of  the  method  to  two  more  events  picking  absolute  minima  of  CUSUM. 
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Combining  CUSUM  with  Simulated  Annealing. 


Instead  of  finding  an  absolute  minimum  in  the  CUSUM-linear  trend  sum,  the  minimum  could 
also  be  located  by  some  randomized  search  method.  Simulated  annealing  is  such  an 
optimization  technique  which  is  designed  to  find  global  minima  of  irregular  functions  where 
many  local  minima  may  exist.  It  tends  to  disregard  minor  local  minima  and  converge  to  the 
lowest  points.  It  uses  the  random  Metropolis  search  algorithm  which  is  based  on  a 
thermodynamic  analogy  ( Press  et  al  1986).  Initially,  it  allows  the  search  using  large  steps  in  the 
independent  variables  which  may  even  be  associated  with  increased  values  of  the  function.  This 
allows  the  solution  to  “jump  out”  from  local  minima  and  resume  search  for  other  minima.  As 
“cooling”  occurs  such  steps  are  accepted  less  and  less  and  finally  the  solution  will  settle  in 
broad  global  minima. 
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Figure  37.  Very  tight  arrival  time  pick  clusters  were  obtained  using  repeated  applications  of  simulated 
annealing.  These  are  processing  results  on  ARCESS  data  for  Kola  events.  Time  length  102.4  seconds. 


Since  there  is  a  possibility  that  occasionally  local  minima  may  be  located,  a  few  rules  are 
needed  to  restart  the  algorithm.  The  algorithm  is  extremely  fast  and  can  be  run  many  times  to 
narrow  down  on  the  actual  arrival  times  and  assess  their  mean  error  empirically  (Figure  37). 
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The  examples  in  Figure  37  are  displays  of  histograms  of  300  independent  SA  searches  vs.  the 
original  and  AR  filtered  seismograms.  We  are  attempting  to  pick  the  emergent  Pn  and  Pg 
phases  on  SPZ  seismograms  from  the  CUSUM.  Note  that  most  clusters  of  the  picks  correspond 
to  where  a  seismologist  would  see  arrivals.  The  variability  in  individual  clusters  of  solutions  can 
be  used  to  assess  the  average  scatter  in  travel  time  picks.  These  are  processing  results  on 
ARCESS  data  for  Kola  events. 
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Figure  38.  Very  tight  arrival  time  pick  clusters  were  obtained  using  repeated  applications  of  simulated 
annealing.  These  are  processing  results  on  ARCESS  data  for  Kola  events.  Time  length  102.4  seconds. 
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Figure  39.  Some  noisy  regional  events  have  multiple,  ambiguous  starts  for  the  main  arrival  groups.  This  is 
reflected  in  the  SA  results  which  form  multiple  clusters.  It  must  be  noted,  however,  that  a  human  analyst 
would  also  have  difficulty  in  picking  onset  times  in  such  cases.  The  Sn  and  Lg  groups  in  this  example 
have  several  energetic  portions  rather  than  a  simple,  impulsive  start.  Time  length  102.4  seconds. 
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Evaluation  of  the  Performance  of  the  CUSUM  Procedure  for  Estimating  Pn  Onset  Times. 


In  the  following  we  evaluate  the  performance  of  the  combined  CUSUM-SA  procedure  for 
picking  onset  times  of  first-arriving  Pn  phases.  The  evaluation  is  based  on  comparing  the 
performance  of  human  analysts  to  a)  picking  the  absolute  minimum  of  CUSUM  formed  around 
the  first  arrival  b)  picking  the  median  of  multiple  picks  using  SA  on  the  CUSUM.  Prefiltering 
preceded  all  the  data  processing.  The  two  kinds  of  prefilters  applied  to  the  seismograms  were  2- 
7  Hz  Butterworth  band-pass  filters  and  filters  designed  by  taking  the  S/N2  spectral  ratio  such 
that  the  maximum  was  set  at  unity  and  cutoffs  were  placed  at  the  values  at  0.24.  The  latter  are 
similar  to  the  filters  that  define  “useful  bandwidth”.  We  have  seen  little  difference  in  the 
performance  of  these  filters  in  accordance  with  the  comments  made  by  Kvaerna  (1996),  The 
events  used  had  originally  very  high  S/N  ratios,  especially  on  the  prefiltered  traces.  To  provide  a 
“true”  onset  time  the  practically  noise-free  original  trace  was  picked  by  the  analyst.  In  order  to 
construct  noisy  data,  we  have  fitted  a  15-th  order  AR  model  to  the  noise  prior  to  the  signal 
arrival  and  this  model  was  used  to  construct  independent  noise  samples  by  filtering  Gaussian 
white  noise.  We  have  verified  that  the  resulting  artificial  noise  samples  had  spectra  very  close  to 
those  of  the  actual  noise  and  thus  the  order  of  AR  process  was  appropriate.  The  S/N  ratio  was 
defined  as  SNR=max(signaI  amplitude)/(Twice  the  std.  deviation  of  noise).  This  value  is 
comparable  to  the  SNR  definition  of  Kvaerna  who  used  a  max(signal)/max(noise)  ratio  with 
short  noise  samples  of  a  few  hundred  points,  and  we  expect  that  the  differences  between  these 
two  definitions  are  inconsequential.  Actually,  our  simulated  records  may  be  noisier  for  the  same 
SNR  value  as  Kvaerna’s,  since  5%  of  the  points  in  a  record  must  exceed  two  standard  deviations 
so  that  the  maximum  in  the  noise  sample  will  be  larger  than  that  value.  The  procedures 
described  here  were  implemented  in  MATLAB  which  provides  a  convenient  environment  for 
quick  prototyping  at  the  expense  of  longer  running  times  of  the  MATLAB  interpreter.  Practical 
software  implementations  of  this  methodology  would  require  the  use  of  the  C  codes  that  we 
have  also  developed  independently,  or  the  compilation  of  the  MATLAB  code  (using  the 
recently  introduced  MATLAB  compiler). 

The  process  goes  through  a  set  of  SNR  values  in  increasing  order  starting  with  the  worst  value 
which  is  0.5  and  increasing  to  9.  In  order  to  prevent  the  analyst  from  locking  into  the  arrival 
time  by  comparing  seismograms  with  varying  SNR  and  by  being  able  to  ‘peek  through’  low 
noise  sections,  we  delayed  the  signal  by  random  amounts  before  presenting  each  simulation. 
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The  simulated  noisy  seismograms  were  presented  to  the  analysts  who  were  asked  to  click  on 
their  best  estimates  of  onset.  Figure  38  shows  some  examples  of  these.  Note  that  the  visible 
spectral  differences  between  signals  and  noise  became  small,  especially  at  low  to  moderate  SNR 
values.  Following  this  the  CUSUM  was  computed  and  the  times  were  picked  automatically. 
Finally,  the  noise-free  trace  is  presented  to  the  analysts  for  their  best  onset  pick  that  was  used  as 
‘true’  onset  time.  The  differences  between  this  value  and  the  other  onset  time  estimates  were 
plotted  for  all  the  three  methods  against  the  logarithms  of  SNR  values,  with  the  random  time 
shifts  corrected  for,  of  course.  This  is  the  same  kind  of  evaluation  method  as  the  one  used  by 
Kvaerna  (1996a,  b)  and  Yokota  et  al  (1981). 

As  the  process  starts  at  low  S/N  there  are,  naturally,  some  complete  failures  in  picking  the 
arrival  time.  We  advised  the  analysts  to  click  close  to  the  ends  of  seismograms  if  they  decided 
that  it  is  not  possible  to  see  the  signal.  This  resulted  in  error  estimates  larger  than  2  sec,  which 
we  defined  as  failures.  The  failures  of  automatic  processing  results  that  followed  were  defined 
similarly  (At  >  2  sec).  The  reader  may  redefine  the  failures  differently  by  setting  this  tolerance 
to  a  lower  level  he  or  she  prefers. 
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Figure  40. 
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Figure  40.  Onset  time  errors  with  simple  3-14  Hz  bandpass  prefiltering. 


The  procedures  were  repeated  with  S/N^  filtering  instead  of  fixed  band  pass  filtering  with  onset 
time  errors  of  similar  accuracy  as  seen  in  the  panels  in  Figure  41.  Apparently,  the  kind  of 
prefiltering,  as  long  as  it  is  reasonable,  does  not  make  much  difference. 
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Figure  4 1 .  Onset  time  errors  the  three  procedures  with  S/N2  filtering. 
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The  surprising  finding  of  these  tests  is  that  most  automatic  onset  times  at  SNR  levels  above  two 
were  comparable  in  quality  to  that  of  the  analyst.  This  is  a  lower  level  than  the  SNR  of  5  quoted 
by  Kvaerna  for  the  current  versions  of  the  AR  processor  at  IDC.  The  reasons  for  this  difference 
are  not  obvious.  It  appears  that  while  the  AR  method  is  sensitive  to  both  RMS  amplitude  and 
spectral  content  changes,  for  the  appropriately  prefiltered  seismograms  the  process  is  primarily 
driven  by  changes  in  the  RMS  amplitude  level.  Spectral  differences  between  the  noise  and 
signal  windows  become  minor  after  such  prefiltering.  Assuming  that  this  is  so,  it  seems  to  be 
computationally  wasteful  to  compare  the  RMS  fits  to  two  not-too-different  AR  models,  i.e. 
spectral  contents,  with  the  associated  instabilities  in  computations,  instead  of  simply  comparing 
changes  in  the  RMS  amplitude  levels  which  is  much  simpler.  The  other  difference  is  that  we 
have  given  up  trying  to  find  the  precise  point  at  which  the  change  occurs  between  the  two  AR 
models.  Rather,  we  try  to  find  regions  where  the  amplitude  increases  occur,  and  accept  the 
ambiguities  associated  with  such  determination.  The  process  will  have  to  be  reevaluated  against 
improved  versions  of  the  AR  processors  at  IDC  currently  under  development. 

Application  of  the  CUSUM-SA  Method  to  Segment  Complete  Regional  Seismograms. 

In  the  previous  examples,  we  applied  the  CUSUM  based  onset  time  picker  to  the  Pn-Pg  and  Sn- 
Lg  arrival  pairs  separately.  Since  the  amplitudes  of  the  seismic  trace  are  much  larger  for  the  last 
two  it  is  difficult  to  subtract  a  single  linear  trend  such  that  all  four  arrivals  coincide  with  minima 
of  the  resulting  function.  Instead  of  subtracting  a  single  linear  trend  one  can  use  other  ways  to 
accomplish  this.  In  our  examples  below  we  have  subtracted  a  two-piece  linear  trend  with  a 
breakpoint  between  Pg  and  Sn.  Other,  equally  effective,  possibilities  include  the  subtraction  of  a 
quadratic  curve  or  some  long-term  average. 

The  basic  assumption  for  the  application  of  such  methodologies  is  that  we  already  have  detected 
the  event  (by  STA/LTA)  and  know  the  approximate  times  of  the  main  arrivals  (from  slowness 
and  polarization  analyses  such  as  those  used  in  IMS  or  the  generalized  versions  of  them  as 
described  by  Der  et  al  1993.  Thus,  we  can  position  the  breakpoint  in  the  trend  between  Pg  and 
Sn. 

As  we  have  indicated  above,  the  various  cluster  analysis  methods  can  be  combined  with 
repeated  applications  of  SA  to  provide  estimates  to  both  the  onset  times  and  their  uncertainties. 
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We  have  applied  the  K-means  clustering  algorithm  to  the  populations  of  independent  onset  time 
estimations  from  SA  ( Tou  and  Gonzalez  1974).  This  is  really  not  the  best  method  to  do  this  task, 
but  since  we  already  had  a  working  program  for  the  K-means  algorithm  it  was  expedient  to  use 
it.  The  1-D  version  of  K-means  was  modified  by  squaring  the  distance  measure,  which  would  be 
otherwise  just  the  linear  distance  between  points.  This  enhanced  the  property  of  the  algorithm  to 
identify  the  tight  clusters  of  onset  times,  instead  of  merging  several  distinct  clusters  into  one. 
When  the  K-means  algorithm  was  used  it  was  found  that  the  best  way  to  apply  this  method  of 
cluster  analysis  was  to  specify  a  larger  number  of  clusters  than  the  number  of  expected  arrivals 
and  discard  the  clusters  with  membership  less  than  10%  of  the  total  onset  time  estimates.  This 
approach  eliminated  the  few  scattered  values  and  concentrated  on  the  few  (typically  four)  main 
arrival  clusters  within  the  search  window. 

Once  the  cluster  centers  are  identified  (K-means  supplies  these),  some  adjustment  for  the  late 
bias  of  the  cluster  center  can  be  accomplished  by  moving  the  estimated  onset  time  by  1.5-2 
standard  deviations  (of  the  cluster)  to  earlier  times.  This  strategy  is  sufficient  for  removing 
visible  bias  due  to  the  shifts  in  minima  caused  by  the  imposition  of  linear  trend  and  thus  put  the 
onset  time  to  where  a  human  analyst  would  pick  it.  Since  onset  times  of  secondary  arrivals  are 
generally  poorly  defined,  the  algorithm  performs  not  worse  than  a  human  analyst.  The  picking 
of  Pn  onset  times  is  a  separate  issue  requiring  a  smaller  time  window. 

Figures  42  to  44  show  examples  of  the  onset  time  determinations  by  the  automatic  algorithm 
described  above  as  applied  to  the  events  analyzed  by  Der  et  al  1993.  In  Figure  42,  we  plot  all 
the  main  functions,  the  raw  trace,  the  4-10  Hz  band-pass  filtered  trace,  the  histograms  of  the  SA 
onset  time  picks  and  the  cumulative  sum  (CUSUM)  of  the  absolute  amplitudes  of  the  band-pass 
filtered  trace.  The  final  onset  time  picks  are  the  vertical  thin  lines.  The  figure  shows  the  results 
for  a  set  of  K2  mine  events.  Since  these  events  have  very  clear  phase  beginnings  and  occurred  at 
the  same  quarry,  the  phases  were  picked  fairly  consistently  from  event  to  event,  i.e.  the  spacings 
of  the  various  automatically  identified  arrivals  are  very  close.  Only  one  Sn  and  one  Pg  phase 
was  missed. 

Similarly  good  results  were  obtained  for  some  other  mines  (Figure  43).  We  have  found  that  in 
many  cases  where  there  were  only  three  main  visible  arrivals  the  algorithm  adjusted  to  this 
situation  automatically. 
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On  the  other  hand,  some  results  were  not  as  good,  especially  some  K.4  mine  events.  The 
waveforms  of  these  events  have  numerous  amplitude  fluctuations  and  the  initial  P  phases  have 
emergent  beginnings.  The  K-means  algorithm  often  groups  scattered  picks  into  arbitrary  groups 
and  the  emergent  beginnings  are  missed.  As  we  have  pointed  out  above  some  of  these  problems 
could  be  eliminated  by  the  substitution  of  some  other  cluster  identification  algorithm  in  place  of 
K-means. 
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Figure  42.  These  are  the  results  for  the  set  of  K2  events  analyzed  by  Der  et  al  (1993).  Since  these  have 
well-defined  phases  with  abrupt  beginnings,  all  were  picked  with  the  exception  of  one  (weak)  Sn  and  a 
Pn.  The  relative  times  of  the  picks  are  also  consistent.  Time  lengths  are  102.4  seconds. 
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Figure  43.  These  are  results  for  three  K5  mine  events  and  two  K8  mine  events.  Time  lengths  are  102.4 
seconds.  Although  some  Pn  times  are  late,  such  errors  can  be  corrected  by  using  separate  procedures  for 
Pn  first  arrivals  as  shown  above. 
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Figure  44.  These  examples  from  the  K4  group  illustrate  cases  where  we  still  had  problems.  These  include 
missing  weak  phases,  missing  emergent  phases,  late  picks  missing  emergent  beginnings  and  having 
double  picks  at  times  where  an  analyst  would  not  pick  anything.  Most  of  these  problems  appear  to  be 
associated  with  the  K-means  grouping  which  will  be  eliminated  in  the  future.  The  main  problems  with 
automatic  onset  time  picking  seemed  to  be  associated  with  the  K-means  algorithm  which  could  be 
eliminated  in  the  future  and  some  more  robust  method  substituted  for  it.  Although  it  seemed  to  work  well 
for  events  sets  with  a  few  well-defined  arrivals,  some  problems  were  encountered  for  events  where  the 
amplitudes  fluctuated,  creating  many  small  clusters.  In  some  cases  sets  of  small  clusters  were  grouped  into 
one  by  the  K-means  algorithm,  thus  creating  an  arrival  that  would  not  be  seen  by  an  analyst  as  such.  We 
see  these  problems  as  due  to  the  K-means  method  and  not  to  the  basic  approach.  The  K-means  algorithm 
will  be  substituted  by  a  cluster  finder  that  progresses  with  time  and  is  combined  with  some  algorithms  that 
limit  the  cluster  size,  spread  and  automatically  eliminates  outliers. 


In  these  tests  we  have  paid  no  attention  to  computational  efficiency;  the  main  goal  was  to 
demonstrate  the  feasibility  of  computer  onset  time  estimation.  The  process  with  300 
individually  estimated  onset  times  and  K-means  sorting,  which  takes  two  seconds  on  a  166MHz 
Pentium  with  compiled  FORT AN  code,  can  be  speeded  up  considerably  by  the  means  described 
next.  First  of  all,  the  starting  points  for  SA  searches  could  be  moved  closer  to  the  known 
approximate  onset  times,  instead  of  being  distributed  randomly  as  above  and  the  initial 
‘temperature’  could  be  lowered  and  fewer  cooling  steps  could  be  taken.  The  cooling  sequence 
could  be  optimized  also.  These  steps  could  easily  reduce  the  amount  of  random  searching  by 
factors  of  two  to  three.  The  added  benefit  for  K-means  would  be  that  the  input  could  be  sorted 
into  fewer  clusters  since  there  would  be  fewer  scattered  onset  times.  The  independent  SA 
searches  can  also  be  parallelized  on  suitable  computers,  thus  reducing  the  time  to  a  fraction  of 
it.  Instead  of  1-D  K-means,  more  efficient  and  simpler  methods  can  be  used  for  identifying  the 
main  clusters.  It  is  certain  that  the  processes  tested  described  above  are  needlessly  computer¬ 
intensive  and  complex  and  thus  could  be  simplified  considerably. 
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Regardless  of  the  complexity  of  the  algorithms  tested,  even  a  two-second  running  time  would  be 
no  obstacle  to  using  the  algorithms  unmodified,  since  the  process  can  always  be  run  in  the 
background,  while  the  analyst  is  busy  with  other  tasks,  and  in  the  fully  automatic  mode  a  two- 
second  running  time  is  sufficiently  fast  for  processing  all  events  in  real  time. 

Multichannel  Deconvolution 

With  regards  to  the  multichannel  deconvolution  method,  we  attempted  to  modify  it  in  order  to 
make  more  useful  for  isolated  teleseismic  events  as  opposed  to  events  clustered  at  test  sites.  The 
main  purpose  is  to  reduce  site-source  tradeoffs  in  analyzing  the  data  ( Der  et  al  1986,  1987, 
1992).  The  latter  scenario  is  more  appropriate  to  isolated  violations  of  a  test  ban  treaty  and  to 
discrimination  of  nuclear  explosions  from  earthquakes  on  a  global  basis.  Applying  the  concepts 
of  our  multichannel  deconvolution  to  such  events  as  seen  at  smaller  arrays  and/or  networks 
requires  some  revision  of  the  procedures  applied.  Since  there  are  no  multiple  events  in  the  same 
source  area,  the  site  spectral  factors  cannot  be  improved  by  event-averaging.  Sites  in  a  network 
with  the  strongest  site  distortion  (non-transparent  sites)  must  be  identified  and  eliminated  when 
estimating  the  source  factors. 

On  the  other  hand,  closely  grouped  sites,  such  as  those  at  small  regional  arrays,  will  have 
similar  site  effects  that  may  influence  the  source  estimates.  Obviously,  we  are  interested  in  the 
sources  and  not  in  the  site  effects.  The  remedy  we  recommend  is  to  test  the  teleseismic  P 
waveforms  for  clustering  using  the  K-means  algorithm.  In  a  large  network  the  cluster  analysis 
will  identify  outlier  stations  (non-transparent  stations).  If,  on  the  other  hand,  the  waveforms 
cluster  according  to  their  geographical  proximity  (such  that  the  waveforms  from  each  regional 
array  would  form  a  cluster),  then  one  should  use  only  a  few  waveforms  from  each  such  array  to 
avoid  contamination  from  site  effects.  We  have  written  a  program  to  perform  K-means  cluster 
analysis. 

Multichannel  deconvolution  for  the  single  event  case  boils  down  to  something  resembling  the 
“phaseless  seismogram”  method  developed  by  Stewart  and  Douglas  1983.  During  the  reporting 
period,  we  have  implemented  the  phaseless  seismogram  as  an  option  that  can  be  used  for  single 
events. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  report  an  adaptive  approach  for  automatic  identification  of  regional  arrivals  is  tested.  The 
work  was  aimed  at  developing  regional  phase  identification  and  detection  techniques  suitable  to 
small  regional  arrays  with  much  fewer  sensors  than  those  present  at  NORESS,  ARCESS  or 
GERESS. 


Figure  45.  A  general  adaptation  scheme  for  recognizing  regional  events. 

The  general  approach  is  to  collect  array  data  for  a  set  of  training  events  (learning  set)  known  to 
be  close  to  each  other  and  compute  a  set  of  digital  filters  that  extrapolate  the  motion  at  one 
sensor  from  that  at  the  other  sensors  for  each  regional  arrival.  The  success  of  this  extrapolation, 
for  all  arrivals  simultaneously,  is  a  criterion  for  the  assignment  of  other  events  to  the  same 
location  (Figure  45).  This  approach  apparently  works  for  the  Kola  peninsula  data  set  we  have 
tested.  The  method  can  be  made  fully  automatic  for  small  arrays  of  the  type  that  will  be  used  in 
global  monitoring.  Despite  these  results,  the  methodology  must  be  tested  using  more  difficult 
cases. 

Extensive  testing  of  the  original  frequency  domain  formulation  presented  by  Der  et  al  (1993) 
was  done.  From  the  work  done  thus  far  we  conclude  that  it  works  better  in  regional  arrival 
discrimination  than  the  conventional  method  based  on  3-component  recti  linearity  and  slowness. 
The  technique  works  well  for  ARCESS  recordings  of  the  K2  and  K4  mines  in  the  Kola  for 
which  9  and  1 1  events  were  available  respectively.  In  these  tests  we  have  used  only  six  sensors, 
a  three  component  combination  at  AO  and  three  vertical  sensors  from  the  C  ring.  We  have  also 
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tested  the  leave-one-out  method  for  both  arrays  for  the  same  combination.  These  tests  showed 
that  the  method  is  very  robust  for  Pn  and  Pg,  but  one  outright  failure  of  picking  Lg  occurred  for 
the  K2  mine  and  some  of  the  correlations  were  significantly  lower  for  Sn  and  Lg.  This  indicates 
that  the  method  will  not  work  well  if  the  number  of  training  events  is  much  below  eight  or  nine. 

For  the  frequency  domain  case  we  have  also  formulated  a  testing  methodology  based  on  the  F 
statistic.  This  method  for  testing  has  not  been  implemented  yet,  as  the  work  in  this  paper  used 
the  normalized  correlation  coefficients  in  Der  etal  (1993). 

Fully  automatic  processing  of  regional  seismograms  appears  to  be  quite  feasible.  After  event 
detection  by  conventional  means  (STA/LTA),  conventional  or  generalized  polarization  analyses 
(Der  et  al  1993)  can  be  used  for  locating  the  main  wave  groups  in  the  regional  seismograms.  In 
the  case  where  no  suites  of  collocated  master  events  are  available,  these  can  be  furnished  by 
standard  slowness-rectilinearity  analyses. 

CUSUM  based  algorithms  can  be  used  effectively  for  determining  the  onset  of  Pn  with  an 
accuracy  comparable  to  a  human  analyst.  Late  bias  for  emergent,  noisy  arrivals  is  unavoidable 
and  common  to  both  the  human  analysts  and  any  automatic  process  and  thus  not  the  sole 
property  of  automatic  algorithms  ( Kvaerna  1995). 

The  onset  times  of  secondary  phases  (Pg,  Sn  and  Lg)  can  be  determined,  somewhat  less 
accurately,  by  CUSUM  based  searches  using  cluster  analysis  on  suites  of  picks  generated  by 
using  repeated  applications  of  simulated  annnealing.  The  algorithms  tested  in  this  report  can  be 
considerably  simplified  in  the  future. 

Full  automatic  analysis  of  regional  seismograms  can  be  accomplished  in  real  time  with 
moderate  amount  of  computer  power  (on  the  486  or  Pentium  PC  level).  S/N  limitations  are 
common  to  all  processes. 

Onset  time  estimates  from  frequency  prefiltering  combined  with  various  forms  of  amplitude 
related  CUSUM  based  algorithms  give  results  that  compare  favorably  with  those  produced  by 
human  analysts  and  the  AR  model  based  onset  time  algorithms. 
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The  exact  nature  of  the  prefilter  for  enhancing  the  SNR  contrast  between  the  arrivals  and  the 
background  noise,  fixed  band-pass,  AR  or  the  more  exact  S/N^,  does  not  seem  to  be  critical  as 
long  as  they  are  reasonably  close  to  the  last  of  these.  Various  parts  of  this  report  used  all  three 
of  these  prefiltering  methods. 

As  revealed  by  examining  regional  signals  with  high  S/N  ratios,  many  of  them  build  up 
gradually  and  have  emergent  beginnings.  Consequently,  late  bias  of  onset  times  for  noisy 
regional  signals  is  common  to  all  automatic  methods  and  the  analyst  determinations  as  well. 

Although  the  method  of  picking  the  minimum  of  the  CUSUM  combined  with  a  linear  trend 
gives  fairly  accurate  onset  times,  it  seems  to  be  more  advantageous  to  utilize  simulated 
annealing  (SA)  for  repetitively  finding  the  minima.  This  approach  gives  populations  of  multiple 
onset  times  that  can  be  further  processed  by  identifying  the  largest  clusters  of  these  as  arrivals 
and  derive  measures  of  their  statistical  uncertainties.  More  work  is  needed  to  define  the  best 
non-parametric  statistical  approaches  to  do  this. 

The  CUSUM  based  methods  are  especially  suited  for  picking  onset  times  of  later  arrivals  in  the 
regional  seismograms.  Since  the  spectral  content  of  a  regional  seismogram  changes  little  as  the 
time  progresses,  onset  time  estimation  must  be  mostly  based  on  amplitude  changes,  not  spectral 
changes  as  in  the  AR  model  based  methods. 

CUSUM-based  methods  to  pick  phase  onset  times  can  also  be  developed  based  on  statistics  that 
are  diagnostic  of  polarization,  slowness,  and  spectral  changes. 
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